#include "Pch.h"

#include "SphericalHarmonics.h"
#include "Common.h"

namespace IS
{

double SphericalHarmonics::P(int L, int M, double x)
{
	double pmm = 1.0;

	if(M>0)
	{
		double somx2 = sqrt((1.0-x)*(1.0+x));
		double fact = 1.0;
		for(int i=1; i<=M; ++i)
		{
			pmm *= -fact * somx2;
			fact += 2.0;
		}
	}

	if(L==M) return pmm;

	double pmm1 = x * (2 * M + 1.0) * pmm;

    if(L==M+1)
        return pmm1;

	double pll = 0;

	for(int ll=M+2; ll<=L; ++ll)
	{
		pll = ( (2*ll-1.0)*x*pmm1 - (ll+M-1.0)*pmm ) / (ll-M);
		pmm = pmm1;
		pmm1 = pll;
	}

	return pll;
}

double SphericalHarmonics::K(int L, int M)
{
    double tmp = ((2*L+1) * Factorial::eval(L-M)) / (4 * Consts::Pi * Factorial::eval(L+M));
	return sqrt(tmp);
}

double SphericalHarmonics::Y(int L, int M, double theta, double phi)
{
    const double sqrt2 = sqrt(2.0); 

    if(M==0) return K(L,0) * P(L, 0, cos(theta));
    else if(M>0) return sqrt2 * K(L, M) * cos(M*phi) * P(L, M, cos(theta)); 
    else return sqrt2 * K(L, -M) * sin(-M * phi) * P(L, -M, cos(theta));
}

}	// namespace IS